1

シミュレートされたべき乗則テール データの一部について、対数対数軸に CCDF グラフをプロットしたいと考えています。以下は、通常の軸に CCDF グラフをプロットする R コードです。リンクのコードを使用しました: ( CCDFグラフをプロットしますか? )

> load("fakedata500.Rda")
> x<-fakedata500
> f<-ecdf(x)
> f
Empirical CDF 
Call: ecdf(x)
 x[1:500] = 0.50174, 0.50307, 0.50383,  ..., 81.674, 140.63
> plot(f)

以下は ECDF グラフです。

ここに画像の説明を入力

> plot(sort(x), 1-f(sort(x)), type="s", lwd=1)

このコマンドは、CCDF グラフを表示します。

ここに画像の説明を入力

ただし、CCDF グラフを両対数軸にプロットして、下の図のようなグラフを作成したいと思います。

ここに画像の説明を入力

Rでそれを行う方法はありますか？

もしそうなら、CCDFグラフでも線形回帰を行うにはどうすればよいですか? 以下のコマンドを使用してみましたが、うまくいきません。

a<-plot(sort(x), 1-f(sort(x)), type="s", lwd=1)
> a
NULL
> res=lm(a)
Error in terms.formula(formula, data = data) : 
  argument is not a valid model

とても感謝しています。

アップデート：

@BondedDust から提供されたコードを使用して、CCDF グラフを正常に生成しました。

(plot(sort(x) , 1-ecdf(x)(sort(x) ), log="xy"))

ここに画像の説明を入力

以下は、データセットを生成する方法のコードです。

u<-runif(500)
fakedata500<-((2*(1-u))^(-1))
4

3 に答える 3

3

分位数を使用する別の方法を次に示します。

library(VGAM)  # for rpareto(...)
set.seed(1)    # for reproducible example
X <- rpareto(1000,location=1,shape=1)
p <- ppoints(100)
par(mfrow=c(1,3))
plot(quantile(X,p=p),p,type="l",ylab="P(X < x)",xlab="x",main="CDF")
plot(quantile(X,p=p),1-p,type="l",ylab="P(X > x)",xlab="x",main="CCDF")
plot(log(quantile(X,p=p)),log(1-p),
     ylab="log[P(X > x)]",xlab="log(x)",main="CCDF: log-log")

そして、これが回帰です。

df  <- data.frame(x=log(1-p),y=log(quantile(X,p=p)))
fit <- lm(y~x,df)
summary(fit)
# ...
# Coefficients:
#              Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.039559   0.007584    5.216 1.02e-06 ***
# x           -0.944380   0.005427 -174.028  < 2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.05317 on 98 degrees of freedom
# Multiple R-squared:  0.9968,  Adjusted R-squared:  0.9967 
# F-statistic: 3.029e+04 on 1 and 98 DF,  p-value: < 2.2e-16
于 2014-05-04T16:50:47.250 に答える
2

[poweRlaw][1]これにはパッケージを使用する必要があります。

パッケージをロードしてデータを生成します。

library(poweRlaw)
x = rplcon(1000, 1, 2)

次に、連続べき乗則オブジェクトを作成します

m = conpl$new(x)

プロット

plot(m)

スケーリング パラメーターを推定します。

p = estimate_pars(m)
m$setPars(p)

適合線をプロットに追加します

lines(m, col=2)

詳細については、パッケージのビネットを参照してください。

ところで、スケーリング パラメーターを推定する際の単純な線形回帰は避けてください (すべての仮定が壊れています)。

于 2014-05-04T17:35:23.000 に答える