23
12月

R语言常用统计分布的蒙特卡洛模拟

为加深对各种常用统计分布的理解,更好掌握R语言对应的各类分布的概率函数(d、p、q、r族)以及广义线性模型的使用,研究了一下常用统计分布数据的模拟生成方法,收获颇多。

各类常用统计分布蒙特卡洛模拟数据生成的大致思路:

1、构造自变量x的均匀分布
2、根据对应分布的均值函数,构造x变量对应的均值。(广义线性模型的link 函数参考
https://en.wikipedia.org/wiki/Generalized_linear_model#Link_function
3、将均值代入,R中对应分布的随机变量生成函数,得到因变量y(例如正态分布为rnorm、泊松分布为rpois)

 

代码:

#norm distribution simulation
set.seed(1234)
num=100
beta0=1
beta1=0.2
x=beta0+beta1*runif(n=num,-1,1)
y=rnorm(num,mean=x,sd=1)
model=glm(y~x,,family=gaussian(link=’identity’))

 

#possion distribution simulation
set.seed(1234)
num=100
beta0=1
beta1=0.2
x=beta0 + beta1*runif(n=num, min=0, max=5)
lambda=exp(x)
y=rpois(n=num, lambda=lambda)
model = glm(y~x, family=poisson(link = log))

 

#Exponential/Gamma distribution simulation
set.seed(1234)
num=100
beta0=1
beta1=0.2
x=beta0 + beta1*runif(n=num, min=0, max=5)
y=rexp(num,rate=exp(-x))
model=glm(y~x,,family=Gamma(link=’log’))
#使用nls模拟
df=data.frame(x,y)
model=nls(y~exp(a+b*x),data=df,start = list(a=0,b=0))

 

#logistic/probit distribution simulation
set.seed(1234)
num=100
beta0=1
beta1=0.2
x=beta0 + beta1*runif(n=num, min=0, max=5)
#logistic distribution logit=log(odds)=log(p/(1-p))
odds=exp(x)
probs=odds/(1+odds)
#probit distribution probit=Cumulative normal pdf
#probs=pnorm(x)

 

y=rbinom(n=num,size=1,prob=probs)
model=glm(y~x1+x2,family = binomial(link=”logit”))

 

#bionimal/Categorical/Multinomial distribution simulation
library(nnet)
y=rbinom(n=num,size=3,prob=probs)
model <- multinom(y ~ x1 + x2)

代码下载


参考资料:《 Monte Carlo Simulation and Resampling Methods for Social Science》

https://www.sagepub.com/sites/default/files/upm-binaries/57233_Chapter_6.pdf

没有评论

第一个在本文留言。

发表评论

名字(必须)
邮箱(不会被公布)(必须)
网址

字体为 粗体 是必填项目,邮箱地址 永远不会 公布。

允许部分 HTML 代码:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>
URLs must be fully qualified (eg: http://www.yeeach.com),and all tags must be properly closed.

超出部分系统将会自动分段及换行。

请保证评论内容是与日志或 Blog 内容相关的,灌水、攻击性或不恰当的评论 可能 会被编辑或删除。

    RSS订阅

    近期文章

    近期评论

    文章归档

    分类目录