2016-10-04 11 views
1

私はヒープ法の定数を推定しようとしています。`nls`は私のモデルのパラメータを見積もることができません

# Function for Heaps law 
heaps <- function(K, n, B){ 
    K*n^B 
} 
heaps(2,117795,.7) #Just to test it works 

のでn = Word Occurrences、およびKBは、個別の私の予測を見つけるために定数でなければなりません値は以下のとおりです。私は次の機能を構築次に

Number of novels DistinctWords WordOccurrences 
1    1   13575   117795 
2    1   34224   947652 
3    1   40353   1146953 
4    1   55392   1661664 
5    1   60656   1968274 

: 私は、次のデータセットnovels_colectionを持っています言葉

私はこれを試してみましたが、それは私にエラーを与える:

fitHeaps <- nls(DistinctWords ~ heaps(K,WordOccurrences,B), 
    data = novels_collection[,2:3], 
    start = list(K = .1, B = .1), trace = T) 

エラー= Error in numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model

任意のアイデアを、私はこのまたは関数をフィットし、Kの値を取得する方法を修正することができる方法でとB

+0

あなたは何を意味するのですか?どのようなログを変換する必要がありますか?とKのために、私はそれが正でなければならないと思うが、私は確信していない –

+1

...それは線形モデルで解くことができるので、非線形モデルの必要はありません。 :) –

答えて

2

y = K * n^Bの両側でログ変換を実行すると、log(y) = log(K) + B * log(n)が得られます。これはlog(y)log(n)の間の線形関係であるため、線形回帰モデルを使ってlog(K)Bを見つけることができます。 minpack.lmで

logy <- log(DistinctWords) 
logn <- log(WordOccurrences) 

fit <- lm(logy ~ logn) 

para <- coef(fit) ## log(K) and B 
para[1] <- exp(para[1]) ## K and B 
1

我々は非線形モデルを当てはめることができますが、私はそれが(Zheyuanで行われたよう)行います対数変換された変数の線形モデルよりも多くの過剰適合になりやすいだろうと思いますが、我々一部の保有データセットの線形/非線形モデルの残差を比較して実験結果を得ることができます。

library(minpack.lm) 
fitHeaps = nlsLM(DistinctWords ~ heaps(K, WordOccurrences, B), 
        data = novels_collection[,2:3], 
        start = list(K = .01, B = .01)) 
coef(fitHeaps) 
#  K   B 
# 5.0452566 0.6472176 

plot(novels_collection$WordOccurrences, novels_collection$DistinctWords, pch=19) 
lines(novels_collection$WordOccurrences, predict(fitHeaps, newdata = novels_collection[,2:3]), col='red') 

enter image description here

関連する問題