W******r
发帖数: 789 | 1 来自主题: ChineseMed版 - 统计的学问 在微观层次上，一种化学药品进入人体后会发生什么化学反应和带来什么生理变化，没
有人能彻底研究清楚。因此，检验药品的安全性和有效性基本上是用统计的方法。中医
之所以不被主流承认，其中很重要的一个原因也是因为中医缺乏统计分析。那么用统计
方法得出的结论有多可靠呢？我们一起来探讨一下。
(1) H0 VS. H1
在假设检验中，零假设H0是处于被保护的位置的，需要很强的证据才能把它推翻。药品
的安全性通常就是对应于两个组没有显著差别的零假设，处于受保护的位置。有一点经
常被人忽略，那就是如果p-value大于阈值，这个时候的结论应该是没有结论，也就是
既不能得出H0不成立的结论，也不能得出H0成立的结论。事实上，如果p-value是0.1左
右，common sense还是应该倾向于认为H0是不成立的，只是没有明显到可以推翻H0的程
度。
举个例子。一个赌场的庄家邀请你和他玩一个游戏。规则很简单：抛硬币。如果硬币正
面朝上，他给你10块钱；如果硬币反面朝上，你给他10块钱。庄家claim说硬币是均匀
的，正面朝上和反面朝上的概率都是0.5。你是一个精明的人。你决定先在旁边观察一
下别人玩的情... 阅读全帖 |
|
W******r
发帖数: 789 | 2 【 以下文字转载自 ChineseMed 讨论区 】
发信人: WtMaster (关仁隐士), 信区: ChineseMed
标 题: 统计的学问
发信站: BBS 未名空间站 (Tue Oct 9 11:25:40 2012, 美东)
在微观层次上，一种化学药品进入人体后会发生什么化学反应和带来什么生理变化，没
有人能彻底研究清楚。因此，检验药品的安全性和有效性基本上是用统计的方法。中医
之所以不被主流承认，其中很重要的一个原因也是因为中医缺乏统计分析。那么用统计
方法得出的结论有多可靠呢？我们一起来探讨一下。
(1) H0 VS. H1
在假设检验中，零假设H0是处于被保护的位置的，需要很强的证据才能把它推翻。药品
的安全性通常就是对应于两个组没有显著差别的零假设，处于受保护的位置。有一点经
常被人忽略，那就是如果p-value大于阈值，这个时候的结论应该是没有结论，也就是
既不能得出H0不成立的结论，也不能得出H0成立的结论。事实上，如果p-value是0.1左
右，common sense还是应该倾向于认为H0是不成立的，只是没有明显到可以推翻H0的程
度。
举个例子。一... 阅读全帖 |
|
W******r
发帖数: 789 | 3 【 以下文字转载自 ChineseMed 讨论区 】
发信人: WtMaster (关仁隐士), 信区: ChineseMed
标 题: 统计的学问
发信站: BBS 未名空间站 (Tue Oct 9 11:25:40 2012, 美东)
在微观层次上，一种化学药品进入人体后会发生什么化学反应和带来什么生理变化，没
有人能彻底研究清楚。因此，检验药品的安全性和有效性基本上是用统计的方法。中医
之所以不被主流承认，其中很重要的一个原因也是因为中医缺乏统计分析。那么用统计
方法得出的结论有多可靠呢？我们一起来探讨一下。
(1) H0 VS. H1
在假设检验中，零假设H0是处于被保护的位置的，需要很强的证据才能把它推翻。药品
的安全性通常就是对应于两个组没有显著差别的零假设，处于受保护的位置。有一点经
常被人忽略，那就是如果p-value大于阈值，这个时候的结论应该是没有结论，也就是
既不能得出H0不成立的结论，也不能得出H0成立的结论。事实上，如果p-value是0.1左
右，common sense还是应该倾向于认为H0是不成立的，只是没有明显到可以推翻H0的程
度。
举个例子。一... 阅读全帖 |
|
b******a
发帖数: 1470 | 4
80
120
##产生数据####
a<-rnorm(10)
b<-rnorm(10)
c<-rnorm(10)
d<-rnorm(10)
##画图， 坐标轴忽略##
boxplot(a,b,c,d,axes=FALSE)
#放置你需要的横轴坐标#
axis(1,1:4,c(120,220,320,420))
#默认纵轴#
axis(2)
#####
box()
#画均值， 画线很麻烦#
points(1,mean(a),col=2)
points(2,mean(b),col=2)
points(3,mean(c),col=2)
points(4,mean(d),col=2)
希望能抛砖引玉 |
|
z**k
发帖数: 378 | 5 set.seed(0)
x1 <- rnorm(1)
x2 <- rnorm(1)
set.seed(0)
y1 <- rnorm(1)
y2 <- rnorm(1)
then x1 = y1, x2 = y2, while x1 and x2 are pseudo-independent. I guess set.
seed is mostly used for debugging |
|
t*****e
发帖数: 364 | 6 sorry I mean using the package tree, rpart seems to be only for regression
tree
below is some sample code, 写的太仓促难免有bug, 不过你可以大概看一下logic
require(tree)
A_train = matrix(rnorm(20000),nc=20)
Label = rnorm(nrow(A_train))
df_train = data.frame(A_train,Label)
A_test = matrix(rnorm(20000),nc=20)
Label_test = rnorm(nrow(A_test))
df_test = data.frame(A_test,Label)
Prob_all = NULL
for (k in 1:100) {
index = sample(1:length(Label),length(Label),replace = T)
indF = sample(NonSelected,floor(length(NonSele... 阅读全帖 |
|
t*****e
发帖数: 364 | 7 sorry I mean using the package tree, rpart seems to be only for regression
tree
below is some sample code, 写的太仓促难免有bug, 不过你可以大概看一下logic
require(tree)
A_train = matrix(rnorm(20000),nc=20)
Label = rnorm(nrow(A_train))
df_train = data.frame(A_train,Label)
A_test = matrix(rnorm(20000),nc=20)
Label_test = rnorm(nrow(A_test))
df_test = data.frame(A_test,Label)
Prob_all = NULL
for (k in 1:100) {
index = sample(1:length(Label),length(Label),replace = T)
indF = sample(NonSelected,floor(length(NonSele... 阅读全帖 |
|
t**i
发帖数: 688 | 8 enj = rnorm(3, mean=3,sd=1)
M1 = matrix(rnorm(40,mean=0,sd=1),ncol=4,byrow=TRUE)
M2 = matrix(rnorm(40,mean=0,sd=1),ncol=4,byrow=TRUE)
require(boot)
simplex(a=enj,A1=M1[,2:4],b1=t(M1[,1]),A2=M2[,2:4],b2=t(M2[,1]),maxi=TRUE)
Linear Programming Results
Call : simplex(a = enj, A1 = M1[, 2:4], b1 = t(M1[, 1]), A2 = M2[, 2:4],
b2 = t(M2[, 1]), maxi = TRUE)
Maximization Problem with Objective Function Coefficients
x1 x2 x3
2.088811 4.114023 2.513823
No feasible solution could be f |
|
m**c
发帖数: 88 | 9 楼上的解释是对的.
但是如果我同时生成几种分布的随机数,比如:
set.seed(1)
x1=rnorm(10) #生成10个正态分布的随机数
y1=rgamma(10) # 生成10个Gamma分布的随机数
............
set.seed(1)
x2=rnorm(10) #生成10个正态分布的随机数
y2=rgamma(10) # 生成10个Gamma分布的随机数
我想 x1 = x 2 应该是没有什么疑问吧.
但是y1 不一定等于y2吧?
也就是说R里面set.seed()函数仅对其后的第一类随机数生成器(本例子中rnorm())起作
用,对其他的随机数生成函数没有作用(rgamma)? |
|
o****o
发帖数: 8077 | 10 就是奔腾mobile.R对应不同的CPU有多个基于ATLAS的RBLAS.DLL
【http://cran.r-project.org/bin/windows/contrib/ATLAS/】。现在提供了给C2D编译的新RBLAS.DLL，更快一些。我把这个新的跟原来用的给老CPU的比较了一下， C2D的比原配在解SVD上快一倍：
原配 （R version 2.8.1）
> x<-matrix(rnorm(1000^2), ncol=1000)
> system.time(y<-svd(x))
用户 系统 流逝
8.42 0.05 8.50
>
###################
P3/P2:
> x<-matrix(rnorm(1000^2), ncol=1000)
> system.time(y<-svd(x))
用户 系统 流逝
6.50 0.11 6.69
>
###################
P4：
> x<-matrix(rnorm(1000^2), ncol=1000)
> system.time(y<-svd(x))
用户 系统 流逝
6.21 0.07 6. |
|
s*r
发帖数: 2757 | 11 不是了，test yourself
x1 <- c(1,1,0,0,0,0)
x2 <- c(0,0,1,1,0,0)
x3 <- c(0,0,0,0,1,1)
y <- c(rnorm(n=2, mean=0), rnorm(n=2,mean=2), rnorm(n=2, mean=4))
x <- cbind(rep(1, times=6), x1,x2,x3))
xxt <- x %*% t(x)
xtx <- t(x) %*% x
e.xtx <- eigen(xtx)
e.xxt <- eigen(xxt)
V <- e.xtx$'vectors'; round(V, digits=2)
U <- e.xxt$'vectors'; round(U, digits=2)
s.v <-e.xtx$'values' # single values
# all these matrix have rank = 3;
d <- sqrt(diag(s.v[c(1,2,3)]))
U3 <- U[,c(1:3)]
V3 <- V[,c(1:3)]
round(U3 %*% t(U3), ... 阅读全帖 |
|
w******a
发帖数: 25 | 12 Here is an R example to impute one or two missing data in each record:
The data will look like
col1 col2 col3
x
x x x
x x
x x
x x x
x
x x x
...
library(Rlab)
alp = 1
K_delta = 2
len_Y1 = 200
#Sample setting:
#Measurment N_
patient Percent
# 1 12
0.18
# 1 2 4
0.05... 阅读全帖 |
|
q**j
发帖数: 10612 | 13 打听一下，glmnet default给前面加一行constant。但是我的X matrix里面已经有了这
个常数了，能否改变glmnet？还是我自己要改变了？
还有
set.seed(1010)
n=1000;p=100
nzc=trunc(p/10)
x=matrix(rnorm(n*p),n,p)
beta=rnorm(nzc)
fx= x[,seq(nzc)] %*% beta
eps=rnorm(n)*5
y=drop(fx+eps)
px=exp(fx)
px=px/(1+px)
ly=rbinom(n=length(px),prob=px,size=1)
set.seed(1011)
cvob1=cv.glmnet(x,y)
以后
coef(cvob1)的结果和
coef(cvob1,s=cvob1$lambda.min)不一样。
我想第二个应该是我们最后想要的东西，但是第一个是什么？
最后，cv.glmnet里面如何能够取出来那个balance L1和L2 penalty的常数呢？
... ...你好歹应该先知道Lasso和Ridge是什么吧？ |
|
w**********y
发帖数: 1691 | 14 # A simple R code to show the best case, if you think it is kinda against
your intuition.
x1 = rnorm(10000)
y = 0.1*x1+sqrt(0.99)*rnorm(10000)
#a=100
#((0.1+a)*0.1+0.99)/(sqrt((0.1+a)^2+0.99))
x2=y+a*x1
cor(x1,y)
cor(x2,y)
summary(lm(y~x1))
summary(lm(y~x2))
summary(lm(y~x1 + x2)) |
|
z******o
发帖数: 1224 | 15 看correlation并不是有效方法，可以举反例：
t = 500:10000
x = 10*t + rnorm(9501,mean=0,sd=10)
y = 10*sqrt(t) + rnorm(9501,mean=0,sd=10)
x和y的correlation ~ 0.98(他们各自不stationary)， 但是他们的trend本质上不同的。
有另一个idea：对每个组，取一个做为y，其他都是x, 如x1, x2,...,xm.
直接regress y on (x1, x2,...,xm)，包含intercept。虽然xi都不stationary但是
least square方法还是valid，这样做的目的是找出一个linear combination，
使得residual=y - beta*X的variance尽可能的小，所以这样得出的residual是最有可
能是stationary的。
然后对residuals做PP/ADF test，得到P value (可能PP更好点，上面的例子，PP test
的simulation result更好点，adf是错的)。
把这些... 阅读全帖 |
|
s*****n
发帖数: 2174 | 16 你这样当然不一样了. 同样大小的noise, 对于两个model的影响不一样.
第一model本身幅度就大, 相对noise的影响就小.
第二个model本身幅度就小, noise的影响相对就大.
公平的比较应该把第二个model里面的noise term减半.
比如
x <- 0:5
result <- matrix(NA, ncol=2, nrow=1000)
for (i in 1:1000){
y1 <- c(200, 190, 180, 170, 160, 150) + rnorm(6) * 5
y2 <- c(200, 195, 190, 185, 180, 175) + rnorm(6) * 2.5
result[i, 1] <- summary(lm(y1~x))$r.squared
result[i, 2] <- summary(lm(y2~x))$r.squared
}
apply(result, 2, mean)
[1] 0.9453929 0.9474898
两者差不多.
严格来说, noise term 不是严格两倍的关系,
只是近似而已. 严 |
|
l*********d
发帖数: 254 | 17 我的问题是：
其实我不知道怎样写下面这个code， 于是就把所有能放的东西方在里面了，然后自己
运行了一下，什么也没得到。希望大家能给与指点，谢谢！
题目如下：
a) Below is the code that I used to simulate the data.
set.seed(9192)
time<-c(3,6,9,12,15,18,21,24,30,36,48,60)
epsilon<-rnorm(n = length(time), mean = 0, sd = 0.014)
beta0<-1
beta1<--0.001
potency<-beta0 + beta1*time + epsilon
set1<-data.frame(time, potency)
set1
#Batch #2 and #3
set.seed(1991)
epsilon<-rnorm(n = length(time), mean = 0, sd = 0.014)
potency<-beta0 + beta1*time + 0.01 + e |
|
q**j
发帖数: 10612 | 18 这两个的区别在哪里？多谢。
> x = rnorm(1000)
> y = 0.5 + x + rnorm(1000)
> z = lm(y~x)
> a = summary(z)
> a$cov.unscaled
> vcov(z) |
|
s*****n
发帖数: 2174 | 19 这个属于稍微有点点高级的R应用了, 光这么说, 估计LZ未必能够搞定.
我给个example吧, 这里应该包含所有需要的东西了.
for (i in 1:10){
var_name <- paste("X_", i, sep = "")
assign(var_name, rnorm(100))
y_name <- paste("Y_", i, sep = "")
assign(y_name, 3 + 5 * get(var_name) + rnorm(100))
formula_name <- as.formula(paste(y_name, " ~ ", var_name, sep = ""))
print(formula_name)
print(lm(formula_name))
} |
|
S******y
发帖数: 1123 | 20 It works for me...
y<-rnorm(100)
x<-rnorm(100)
fm<-lm(y~x)
plot(y~x)
abline(fm, col=4) |
|
m**c
发帖数: 88 | 21 关于这个set.seed()函数的作用，有点迷糊，哪位给解答一下。
这个set.seed(123456)调用一次，是不是仅仅影响其随后的随机数的生成？
比如：
set.seed()
.....
rnorm(N) # 1
.....
rnorm(N) # 2
是不是set.seed()仅仅影响第一个随机数生成，与第二个无关？
是不是这么个情况？ |
|
w******a
发帖数: 25 | 22 Here is an R example to impute one missing data in each record,half of the code is to make data sample, you probably only need second half,but including them here helps you understand what is going on:
The data will look like
col1 col2
x
x x
x
x x
x x
...
library(Rlab)
alp = 1
Prob_R1 = 0.5
Prob_R0 = 1 - Prob_R1
len_Y1 = 200
K_delta = 2
Y1 = rnorm(len_Y1,mean=0,sd=1)
R1 = rbinom(n=len_Y1, size=1, prob=Prob_R1)
Y2 = rnorm(n=len_Y1,... 阅读全帖 |
|
b*****n
发帖数: 685 | 23 我数学不是很好，老是靠经验，试了试：
eigen(cov(matrix(rnorm(10), nrow = 2)))
eigen(cov(matrix(rnorm(10), ncol = 2)))
您数学牛，给解释一二？ |
|
x****n
发帖数: 43 | 24 sample1 <- rnorm(n1, mean=mean1, sd=sd1)
sample2 <- rnorm(n2, mean=mean2, sd=sd2)
t.test(sample1, sample2) |
|
d******e
发帖数: 7844 | 25 > y = rnorm(100)
> x1 = rnorm(100)
> cor(x1,y)
[1] -1.254646e-05
> x2 = x1 + 0.0000001*y
> cor(x2,y)
[1] -1.243022e-05
> cor(x2-x1,y)
[1] 1
自己试去。 |
|
o****o
发帖数: 8077 | 26 no problem.
you do save a lot of time
> m<-matrix(rnorm(1000^2), ncol=1000)
> ident<-diag(rep(1, 1000))
> a<-rnorm(50)
>
> system.time(ia<-ident%x%a)
用户 系统 流逝
0.95 0.14 1.09
> rm(ia); gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 122044 3.3 350000 9.4 350000 9.4
Vcells 2079203 15.9 85942660 655.7 102081249 778.9
> system.time(res<-m%*%(ident%x%t(a)))
用户 系统 流逝
12.40 0.40 12.85
> rm(res);gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 122048 ... 阅读全帖 |
|
o****o
发帖数: 8077 | 27 pearson correlation is very data specific and can trap you into faked beliefs.
****************
try this code:
set.seed(1)
x<-1:100
y<-0.2*x+rnorm(length(x))*sqrt(var(x))
cor(x,y)
plot(x,y)
x2<-1:100
y2<-900+1.2*x+rnorm(length(x))*sqrt(var(x))
cor(x2,y2)
zx<-c(x,x2)
zy<-c(y,y2) |
|
s*****n
发帖数: 2174 | 28 也不完全是. R在一定程度上是pass by reference/shallow copy的, 从下面这个例子
可以看出来.
> rm(list = ls()); gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 114487 6.2 350000 18.7 350000 18.7
Vcells 122007 1.0 2818617 21.6 3122026 23.9
> a <- matrix(rnorm(1000000), 1000, 1000); gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 114506 6.2 350000 18.7 350000 18.7
Vcells 1122323 8.6 3523602 26.9 3122340 23.9
> b1 <- b2 <- b3 <- b4 <- a; gc()
used (Mb) gc trigger (Mb... 阅读全帖 |
|
s*********e
发帖数: 1051 | 29 > n <- 1000000
> set.seed(2013)
> ldf <- data.frame(id1 = sample(n, n), id2 = sample(n / 100, n, replace =
TRUE), x1 = rnorm(n), x2 = runif(n))
> rdf <- data.frame(id1 = sample(n, n), id2 = sample(n / 100, n, replace =
TRUE), y1 = rnorm(n), y2 = runif(n))
>
> # METHOD 1: MERGE
> system.time(join1 <- merge(ldf, rdf, by = c("id1", "id2")))
user system elapsed
54.028 11.229 65.673
>
> # METHOD 2: PLYR
> # library(plyr)
> # system.time(join2 <- plyr::join(ldf, rdf, by = c("id1", "id2"), type =... 阅读全帖 |
|
a******e
发帖数: 119 | 30 for (g in 1:g)
{
for (m in 2:m)
{
lamda[m,]=sapply(Y[g,1:n],update.lamda)
beta[m]=rgamma(n=1,shape=n*alpha[m-1]+beta.shape,rate = (sum(lamda[m-1,]
)+beta.rate))
temp = rnorm(n=1,mean=alpha[m-1],sd=sigma)
den = (beta[m-1]^(n*alpha[m-1]))*((prod(lamda[m-1,]))^(alpha[m-1]-1))*(
exp(-alpha[m-1])*alpha.rate)/(gamma(alpha[m-1]))^n
num = (beta[m-1]^(n*temp))*((prod(lamda[m-1,]))^(temp-1))*(exp(-temp)*
alpha.rate)/(gamma(temp))^n
accep.prob=num/den
if((accep.prob>=u[... 阅读全帖 |
|
v*******e
发帖数: 11604 | 31
这太容易了。Pvalue can be as small as you want.
N=c(100,1000,10000,100000)
PValue=rep(as.numeric(NA),length(N))
RSquared=rep(as.numeric(NA),length(N))
for (i in 1:length(N)){
sigmma=rnorm(N[i],0,1)
x=rnorm(N[i],0,1)
y=0.1*x+sigmma
RSquared[i]=summary(lm(y~0+x))$"r.squared"
PValue[i]=summary(lm(y~0+x))$"coefficients"[1,4]
}
data.frame(N,RSqaured,PValue)
N RSquared PValue
1 1e+02 0.007324693 3.947863e-01
2 1e+03 0.014429159 1.392340e-04
3 1e+04 0.010822982 1.792036e-25
4 1e... 阅读全帖 |
|
k*******a
发帖数: 772 | 32 多快叫快？
a = matrix(rnorm(2000*2000), nrow=2000)
b = matrix(rnorm(2000*2000), nrow=2000)
system.time(a %*% b)
user system elapsed
4.705 0.004 4.708
如果你的时间远远超过这个，那估计是硬件的问题 |
|
o****o
发帖数: 8077 | 33 你的也太慢了，用的default的BLAS吧
> a = matrix(rnorm(2000*2000), nrow=2000)
>
> b = matrix(rnorm(2000*2000), nrow=2000)
>
> system.time(a %*% b)
user system elapsed
1.47 0.04 0.26 |
|
o****o
发帖数: 8077 | 34 你的啥机器啊，用的R自带的BLAS吧
我的i5-3570k，没超频
上一个是至强的六核CPU
> a = matrix(rnorm(2000*2000), nrow=2000)
> b = matrix(rnorm(2000*2000), nrow=2000)
> system.time(a %*% b)
user system elapsed
1.42 0.14 0.43
> |
|
l****i
发帖数: 398 | 35 如果嫌麻烦自己装open blas什么的, 可以直接装一个revolution r open. RRO自带
Intel Math Kernel Library (MKL)可以multi-thread矩阵运算。我用的i5-3320m, 16g
内存，只能2个线程，速度也还行。
> if(require(RevoUtilsMath)){
+ setMKLthreads(2)
+ }
Loading required package: RevoUtilsMath
>
> a = matrix(rnorm(2000*2000), nrow=2000)
> b = matrix(rnorm(2000*2000), nrow=2000)
>
> system.time(a %*% b)
user system elapsed
0.92 0.00 0.52 |
|
g*r
发帖数: 67 | 36 R的merge效率太低了。根据一个列合并两个dataframe，长的有1m行，短的1k行：
long =data.frame(a=seq(1,1000000), b=rnorm(10) )
short=data.frame(a=seq(1,2000, 2), c=rnorm(10) )
system.time(merge(x=long, y=short, all.x=T))
user system elapsed
10.286 0.370 10.624
4m行就要将近1分钟。8m行超过2分钟。16m要4分钟。
我的真实数据比例子里的复杂些，且有100m行或更高，卡了半天没动静，只好强行终止
进程。对这个效率有点失望。早听说merge很慢，这次领教了。
如果不依靠外部程序（shell，python之类的script），大家是怎么多快好省做合并的
呢？菜鸟一个，轻拍。 |
|
B*O
发帖数: 105 | 37 用match比较快。
> long =data.frame(a=seq(1,1000000), b=rnorm(10) )
> short=data.frame(a=seq(1,2000, 2), c=rnorm(10) )
> system.time(merge(x=long, y=short, all.x=T))
user system elapsed
7.253 0.071 7.385
> system.time(long$c<-rep(as.numeric(NA),dim(long)[1]))
user system elapsed
0.002 0.000 0.001
> system.time(long$c[match(short$a,long$a)]<-short$c)
user system elapsed
0.119 0.000 0.119 |
|
i*r
发帖数: 83 | 38 不想复杂， 你就是要看 10 个feautre之间的关系， 做一个correlation/covariance
matrix， 这个matrix就capture了所有10个feature两两之间相关性。 然后在这个cor
matrix做个hiearchical clustering， 你就看出来相互之间的关系， 这样最简单，
因为correlation matrix 是normalized cov， 阈值范围 （-1,1）， 一幕了然。
R code for simulation variable 7，8,9,10 和variable 1 是正相关：
x = matrix(rnorm(40,mean=2),ncol=10)
for(i in 7:10) x[,i] = x[,1]*i + rnorm(4)
x.cor = cor(x)
gplots::heatmap.2(x.cor)
你会看见1,7-10 会group在一起的， 还不明白的话google “hierachical clustering
” 和 “heatmap”
复杂点的， 做 multidimension scal... 阅读全帖 |
|
d********m
发帖数: 3662 | 39 vec <- c()
for(i in 1:10000){
a <- t.test(matrix(sample(rnorm(10000), 200, rep = FALSE), 100, 2))
vec[i] <- a$p.value
}
有没有异议？ |
|
m********0
发帖数: 2717 | 40 The thing is the it's not a graph for any stock. I just make it with the
following R code,
require(TTR)
ret<-cumsum(rnorm(252,sd=0.036))
ret2<-exp(ret)
p<-50*ret2
ma20<-SMA(p,20)
ema20<-emaTA(p,20)
ema50<-emaTA(p,50)
plot(p, type="l")
lines(ema20, type="l",col='green')
lines(ema50, type="l",col='red')
and I used lognormal distribution for the returns, textbook assumptions
for
a big group of models.
skeptic(as exists everywhere), yes, I borrowed this idea from ET. Please
don
't rush to criticize,... 阅读全帖 |
|
s*******e
发帖数: 664 | 41 ☆─────────────────────────────────────☆
stallone (stallone) 于 (Fri Dec 18 22:24:21 2009, 美东) 提到:
最近有个想法，就是想用C/C++做统计计算，因为用R，matlab都很慢，效率不高。
下面具体说下我的疑惑吧。
我的C的就是学过谭浩强那本C语言，然后就在没有接触其他的了。
而统计计算主要涉及向量，矩阵，三维数组等的运算，当然还有其他的优化问题。
我总是想有些在R里面很容易实现的东西在C里面好像很困难。
如 生存伪正态随机数，在R里就 rnorm 可以实现，
如 矩阵求逆，R 用 solve 即可，还有很多例子。而这些要
用C写的话，就很慢，出错的可能也增多，即使运行起来快。
所以请板上用C/C++做统计计算的人，指点一下，应该怎么入门
用C做统计计算，要注意些什么。
如能推荐相关资料那是再好不过的了。
在此先谢过！
☆─────────────────────────────────────☆
pptwo (pp) 于 (Fri Dec 18 22:45:04 2009, 美 |
|
z**k
发帖数: 378 | 42 我觉得liujz606的答案是正确的，variance应该是T^3/3，你的计算我还没看，等下回
家研究下
#!/usr/bin/R
record <- numeric(10000)
for (i in 1:10000) {
grid <- cumsum(rnorm(1000, sd=sqrt(0.001)))
record[i] <- sum(grid)/1000
}
var(record)
答案是
0.3360694 |
|
b*********n
发帖数: 2284 | 43 你这10*t + rnorm(9501,mean=0,sd=10)是return sequence?
的。 |
|
K*****Y
发帖数: 629 | 44 题目是已经给定可以生成standard normal random variable x, 所以不清楚为什么楼
上的解法还要用到 1/(sqrt(2*pi))和对数变换，感觉用不着这么麻烦。
用R测试并写的程序：
N=10000 # number of random samples, which can be specified
x=rnorm(N) # generate N standard normal random variables
alpha=10 # calculate P(x>alpha)
est=mean(as.numeric(x>0)*exp(-alpha*x-0.5*alpha^2)) #est gives the required
value
思路很简单。要求P(x>alpha)，E[1_{x>alpha}],做一个变换x'=x-alpha，也就是要求
期望值E{1_{x'>0}*w(x')}，其中w(x')=exp(-alpha*x'-0.5*alpha^2) is the weight
function, 也就是新旧变量的likelihood... 阅读全帖 |
|
w**********y
发帖数: 1691 | 45 setwd("C:")
(data <- data.frame(var1=rep(1:2,4), var2=sample(1:100,8), value = rnorm(8)))
(data <- data[with(data,order(var1,var2)),])
level.var1 <- unique(data$var1)
for (v in level.var1){
subdata <- data[data$var1==v,]
write.csv(subdata, paste0("format2",v,".csv"), row.names=F)
write.table(matrix(subdata$value,,2), paste0("format1",v,".csv"), row.
names=F, col.names=F, sep=',')
}
Baozi pls |
|
a******6
发帖数: 78 | 46 一个简单的程序，不知道哪里出了问题，请各位大虾指点一下：
谢谢
#include
#include
#include
#include
int main()
{
double a[3];
int i;
GetRNGstate();
for(i=0;i<=3;i++)
{
a[i]=rnorm(0.0,1.0);
printf("a[%2d]=%f \n",i,a[i]);
}
PutRNGstate();
}
错误如下：
Linking...
test.obj : error LNK2001: unresolved external symbol _PutRNGstate
test.obj : error LNK2001: unresolved external symbol _Rf_rnorm
test.obj : error LNK2001: unresolved external symbol _GetRNGsta |
|
s*****n
发帖数: 2174 | 47 x <- array(1:12, dim=c(2,2,3))
for (i in 1:10){
temp.name <- paste("data", i, sep = "")
temp.data <- data.frame(ss = matrix(rnorm(3)))
names(temp.data) <- paste("name", i, sep = "")
assign(temp.name, temp.data)
} |
|
d*******1
发帖数: 854 | 48 you need to change your mindset when you switch from R, matlab to SAS, which
is a basically data language.
If your variable x is stored in a data set in which variable x has 100
observation. such as
data test;
do i=1 to 100;
x=rnorm(1); output;
end;
run;
One way to call them is to make every observation a macro variable as
follows:
data _null_;
set test;
call symput('x_'||trim(left(put(_n_,best.))),trim(left(put(x,best.))));
run;
Now, you can call any x by using macro |
|
m********h
发帖数: 60 | 49 正在用 R 做simulation,用到 rnorm，所以每次结果都不一样，我就用set.seed（）确
定结果，可是结果很依赖于seed。想问一下是在不同的seed里挑个结果最好的吗？还是
有别的办法？多谢了。 |
|
a********a
发帖数: 346 | 50 I have a program as following,
beta1=2
beta2=8
#for (i in 1:2){
obs=2
x1=matrix(NA,obs,1)
x2=matrix(NA,obs,1)
for (g in 1:obs){
x1[g,]=runif(1,1,2)
x2[g,]=rnorm(1,1)
}
data=cbind(x1,x2)
data
#d=rbind(data[i])
#}
I run the program 2 times, and get the data like following each time,
> data
[,1] [,2]
[1,] 1.452696 0.2718456
[2,] 1.514587 1.2971293
> data
[,1] [,2]
[1,] 1.172726 -0.9674824
[2,] 1.159079 0.5483838
how can I combine these data by row, i.e. I want to get d |
|
