出家如初,成佛有余

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

Posted in 技术相关 by chuanliang on 2017/12/23

为加深对各种常用统计分布的理解,更好掌握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

Advertisements

R中fitted()与predict()函数关系以及predict type参数解析

Posted in 技术相关 by chuanliang on 2017/12/17

R文档对fitted()、predict()函数以及predict函数type参数取值的具体含义说得很不清楚,网上也没有清晰的解释。总结一下。

一、说明

目的一、R中fitted和predict的关系
目的二、以logistic为例,解析predict中type参数不同取值的关系

1、R中fitted和predict的关系
fitted对无link函数模型(例如线性回归),fitted值和predict的值相同
对有link函数的模型(例如logistic、Binomial),
fitted是link函数作用前(或者link逆函数inverse of the link function作用后)的预测值
predict是link函数作用后(或link逆函数inverse of the link function作用前)的预测值
link function  https://en.wikipedia.org/wiki/Generalized_linear_modelLink_function

2、logistic中predict中type参数不同取值的关系
predict中type参数的可选项有:default, type="link",type="response",type="terms"

假设log(odds)=log(p/(1-p))=alpha+beta1*x1+beta2*x2+…+betan*xn
a.、type="link",type为缺省值
  type="link"为缺省值,给出logit线性函数预测值,link=alpha+beta1*x1+beta2*x2+…+betan*xn

b、 type="response"
  type="response" 给出概率预测值。
  response=predict(model,type = "response"),则
  log(response/(1-response))=alpha+beta1*x1+beta2*x2+…+betan*xn=link,或:link=exp(response)/(1+response)

c、type="terms"
  type="terms"表示各个变量的预测值(包括前面参数在模型中的值)。
  terms=predict(model, newdata = traindata,type = "terms")的结果为各个变量term的值
  term1=beta1*x1,term2=beta2*x2, …
  sum(term1,…,termn)+attr(terms,"constant")=log(odds)=log(response/(1-response))
  attr(terms,"constant")的含义 http://blog.minitab.com/blog/adventures-in-statistics-2/regression-analysis-how-to-interpret-the-constant-y-intercept

二、代码

代码下载

#fitted和predict的关系演示
#linear regression
#no link function
num=100
x1=rnorm(num)
x2=rnorm(num)
lr=1+2*x1+3*x2
model=lm(y~x+x2)

fit=fitted(model)
pred=predict(model)
print(all.equal(fit, pred))
#对无link函数的模型,fitted值与predict值相同

#possion
#link=log
num=100
x=rnorm(num)
y=rpois(num,lambda = exp(x))
model = glm(y~x, family="poisson")
fit=fitted(model)
pred=predict(model)
print(all.equal(log(fit),pred))
print(all.equal(fit,exp(pred)))
#对有link函数的模型,fitted值为link函数作用前的值,predict为link作用后的值

#logistic
#link=logit=log(p/(1-p))
num=100
x1=rnorm(num)
x2=rnorm(num)
lo=1+2*x1+3*x2
odds=exp(lo)
probs=odds/(1+odds)
y=rbinom(n=num,size=1,prob=probs)
model=glm(y~x1+x2,family = binomial(link="logit"))

fit=fitted(model)
pred=predict(model)
print(all.equal(log(fit/(1-fit)),pred))
print(all.equal(fit,exp(pred)/(1+exp(pred))))
#对有link函数的模型,fitted值为link函数作用前的值,predict为link作用后的值

#以logistic模型为例子,验证predict中type为缺省值,type="link",type="response",type="terms"关系
#predict中type="link"时的fitted及predict
link=predict(model,type="link")
print(link)
print(all.equal(pred,link))
#type=link为缺省值

#predict中type="response"时的fitted及predict
response=predict(model,type="response")
print(response)
print(all.equal(fit,response))
print(all.equal(pred,log(response/(1-response))))

#type="terms"时的fitted及predict
terms<- predict(model,type="terms")
term=terms[,1]+terms[,2]+attr(terms,"constant")
print(term)
print(all.equal(link,term))
print(all.equal(response,exp(term)/(1+exp(term))))

参考资料:

https://stackoverflow.com/questions/12201439/is-there-a-difference-between-the-r-functions-fitted-and-predict