Original link:tecdat.cn/?p=6534

Original source:Tuo End number according to the tribe public number

 

data

This is a very simplified example. I simulated 1,000 counting observations, with an average of 1.3. Then, if only two or higher observations are observed, I compare the original distribution with the one I got.

 

This code generates:

 

# Raw data:
set.seed(321)
a <- rpois(1000.1.3)

# Truncated version of data:
b <- a[ a > 1]

# graphics:
data_frame(value = c(a, b),
  ggplot(aes(x = value)) +
   (binwidth = 1, colour = "W # model fits the original model well: mean(a) (a,"Poisson") # fitdistr(b,"Poisson")
Copy the code

Estimating the key parameters of lambda complete data (a) worked well, with an estimated value of 1.347, just over one standard error of the true value of 1.3.

Maximum likelihood

Use truncated versions of the dPOis and pPOis functions in Fitdist.

#------------- uses MLE fit ------------------- in R
dtruncated_poisson <- function(x, lambda) {
}
ptruncated_poisson <- function(x, lambda) {
}

fitdist(b, "truncated_poisson", start = c(lambda = 0.5)) 
Copy the code

Note that to do this, I specify a lower threshold of 1.5; Since all data are integers, this actually means that we only observe 2 or more observations. We also need to specify a reasonable starting value for the estimate, lambda, so that the error is not too great.

 

The bayesian

Stan can easily describe data and probability distributions as truncated for the alternative Bayesian approach. In addition to the raw data that I X called in this program, we need to tell it how many observations (n), lower_limit truncation, and any variables needed to characterize the prior distribution of our estimated parameters.

The key parts of the following procedure are:

  • indataOf the specified dataxLower bound for thelower_limit
  • inmodel, specifyxBy truncated distributionT[lower_limit, ]
data {
  int n;
  int lower_limit;
  int  x[n];
  real lambda_start_mu;
  real lambda_start_sigma;
}

parameters {
  reallambda;
}

model {
  lambda ~ normal(lambda_start_mu, lambda_start_sigma);
  
  for(i in 1:n){
    x[i] ~ poisson(lambda) T[lower_limit, ]; }}Copy the code

Here is how R provides data to Stan:

#------------- Call Stan-------------- from R
data <- list(
  x = b,
  lower_limit = 2,
  n = length(),
  lambda_start_sigma = 1
)

fit <- stan("0120-trunc.stan", data = data, cores = 4)


plot(fit) + 
  labs(y = "Estimated parameters") +
  theme_minimal(base_family = "myfont")
Copy the code

Results provide a posterior distribution estimated by lambda and Fitdistrplus methods: 1.35 with a standard deviation of 0.08. Image of confidence interval:

 

Most welcome insight

1.R language Poisson regression model analysis case

2. Numerical simulation in R language: Poisson regression model was simulated

3. Poisson regression analysis for R language

4.R language simulation and dynamic visualization of Buffon needle drop (Buffon needle drop) experiment

5. Use R language to simulate the mixed queuing random service queuing system

6. VaR comparison of GARCH (1,1), MA and historical simulation method

7.R language to do geometric Brownian motion simulation of complex financial products

8. Numerical simulation in R language: Poisson regression model is simulated

9. Pricing of reinsurance contracts under catastrophe risk under R language: Generalized linear Models and Pareto distributions