W niniejszym paragrafie przyjrzymy się kilku popularnym w literaturze mo-delom procesu oceniania produktów oraz momo-delom procesu produk-cji/oceniania. Mają one zastosowanie w Problemie Oceny Producentów.
Będą to dwa modele deterministyczne: model o stałej produktywności i sta-łych przyrostach ocen (rozdz. 4.1.1) oraz model oparty na prawie potę-gowym Lotki (rozdz. 4.1.2), a także niedeterministyczne: model Burrella (rozdz. 4.1.3) oraz, będący w centrum naszych zainteresowań, model ocen i.i.d.(rozdz. 4.1.4).
4.1.1 Model deterministyczny o stałej produktywności i przyrostach ocen Najprostszy, deterministyczny model zachowania się funkcji wpływu dla da-nych o określonej postaci rozważał już sam Hirsch na przykładzie zapropo-nowanego przez siebie indeksu h [107]. W modelu tym zakłada się, że pro-ducent wytwarza stałą liczbę produktów p ∈ N na jednostkę czasu (np. na
rok). Dodatkowo, ocena każdego produktu wzrasta o stałą c 1 po upływie umownego okresu.
Tym samym rozpatrujemy następujący ciąg wektorów x(i) w kolejnych jednostkach czasu i = 1, 2, . . . :
x(1) = (p∗ 0) x(2) = (p∗ c, p ∗ 0) x(3) = (p∗ 2c, p ∗ c, p ∗ 0)
... ... ...
Łatwo zauważyć, że Prod(x(i)) = ip oraz że Sum(x(i)) = pci(i− 1)/2.
Można pokazać [zob. 107], że
H(x(i))≃ pc
p + c(i− 1), (4.1)
zatem wartość indeksu h w tym modelu wzrasta liniowo wraz z czasem.
Ponadto, w pracy [200] uwzględniono (w kontekście bibliometrycznym) analizę wpływu samocytowań na wartość tego indeksu. Warto zwrócić uwagę na fakt, iż jeśli autor cytuje w każdej pracy wszystkie swoje dotychczas opu-blikowane dzieła, nie otrzymując żadnych dodatkowych cytowań, to przy za-łożeniu rocznej produktywności p zachodzi
H(x(i))≃ p
2(i− 1). (4.2)
Modyfikacja tego modelu była też rozpatrywana w [101].
4.1.2 Model deterministyczny oparty na prawie Lotki
W 1926 r. zostało sformułowane tzw. prawo Alfreda Lotki [por. 60, 148],
Prawo Lotki
będące jedną z bardziej znanych zasad (podkreślmy, empirycznych) biblio-metrii. Rozważmy pewną ustaloną grupę twórców. Prawo potęgowe Lotki (ang. Lotka’s power law) głosi, że liczba p(i) autorów, którzy opublikowali i prac, może być opisana funkcją potęgową:
p(i) = K
iα, i = 1, 2, . . . , (4.3) gdzie K i α są dodatnimi stałymi. Najczęściej przyjmuje się α ≃ 2, ale wartość tego parametru może być zależna od konkretnej dziedziny. Prawo to może również służyć np. do opisu liczby cytowań kolejnych uszeregowanych artykułów.
Niech x ∈ In dla pewnego n. Istnieje wzajemnie jednoznaczna odpowied-niość pomiędzy analizowanymi przez nas ciągami ocen produktów a funkcją
p. Można bowiem przyjąć x(n−i+1) = p(i), bowiem p(i) oznacza wartość oceny produktu o i-tej randze. Na podobnej idei bazowała wprowadzona przez nas funkcja ocen produktów πx (por. s. 55).
W teorii tzw. Lotkaian informetrics, która jest od kilkunastu lat konse- IPP kwentnie rozwijana przez L. Egghego, rozpatruje się dla ułatwienia obliczeń uogólnienie funkcji p. Zakłada się wtedy, że p jest ciągłą funkcją określoną na pewnym przedziale liczb rzeczywistych [por. 60]. Charakteryzuje ona wtedy pewien możliwy proces wytwarzania informacji (IPP, ang. information pro-duction process, zob. np. [71]). Jest to podejście w sensie ideowym podobne do naszego, polega ono bowiem na abstrahowaniu od bibliometrycznego kontek-stu przetwarzanych danych. Niemniej, jak widzimy, posługuje się ono innymi metodami opisu.
Jednakże tego typu uogólnienie funkcji p pociąga za sobą problemy in-terpretacyjne. Ponadto, jak widzimy, model ten jest stricte deterministyczny.
Nie da się więc stosować w nim metod wnioskowania statystycznego. Z tych względów podejście to nie budzi raczej większego zainteresowania wśród ba-daczy. Rozpatrywane jest ono głównie przez Egghego i kilku jego najbliższych współpracowników.
I tak w literaturze rozważano m.in. wyznaczenie wartości indeksu h [71]
czy indeksu g [62] dla danych zgodnych z prawem potęgowym, wpływ do-dawania bądź usuwania elementów na wartość h [70], badanie wpływu za-stosowania pewnej klasy transformacji wobec funkcji p na wartość indeksu h [65] czy też wpływu zmiennej czasowej w modelu Lotki na wartości indeksów bibliometrycznych [63, 69].
Przytoczmy teraz opis dwóch ciekawych modeli nieterministycznych.
4.1.3 Model Burrella
W pracach [33, 37–39] Q. Burrell zaproponował następujący stochastyczny model procesu produkcji/oceniania. Zakłada się tu, że oceny poszczególnych produktów są liczbami ze zbioru N0.
Niech t = 0 oznacza chwilę rozpoczęcia działalności producenta.
1. Proces wytwarzania produktów ma rozkład Poissona z parametrem θt, tzn. do chwili t > 0 liczba produktów Yt ∼ Poi(θt). Zauważmy, że z tego założenia wynika, iż czasy pomiędzy wytworzeniem kolejnych produktów można opisać rozkładem wykładniczym o wartości oczeki-wanej 1/θ.
2. Jeśli ti 0 oznacza chwilę powstania i-tego produktu, to ocena roz-patrywanego produktu do chwili t > ti jest zmienną losową z rozkładu Poissona o parametrze Λi(t− ti), tzn. Xi,t∼ Poi(Λi(t− ti)).
3. Intensywność oceny i-tego produktu Λi jest zmienną losową z rozkładu Γ(ν, α), opisanego gęstością
fΛi(x) = αν
Γ(ν)xν−1e−αx.
Model Burrella jest zatem opisany trzema parametrami (por. rys. 4.1). War-tość θ > 0 oznacza oczekiwaną liczbę wytwarzanych produktów na jednostkę czasu, a ν 1, α > 0 określają oczekiwaną intensywność oceny E Λi = ν/α.
Rysunek 4.1: Graficzna reprezentacja modelu Burrella.
Burrell uzyskał kilka ciekawych, acz elementarnych wyników teoretycz-nych, m.in. rozkład wartości łącznej oceny produktów, czyli funkcji Sum [33]
czy też wzór (wyznaczalny numerycznie) na wartość oczekiwaną indeksu h [37], w przybliżeniu zależną logarytmicznie od θ. Badał także zachowanie się innych indeksów bibliometrycznych, np. g i hα [40]. Jednakże, model ten nie doczekał się jeszcze weryfikacji empirycznej poza analizą przebiegu kariery jednego naukowca [39]. Mimo tego, model Burrella jest interesujący, ponieważ stanowi próbę uwzględnienia zarówno zmienności wartości ocen produktów, jak i liczby tych produktów.
W niniejszej dysertacji jednak skupimy uwagę na szczegółowej analizie podstawowych własności probabilistycznych i statystycznych wybranych ope-ratorów agregacji przy założeniu z góry ustalonej liczby produktów. Tego typu analiza nie została do tej pory przeprowadzona w literaturze. Ewen-tualne rozszerzenie tego zagadnienia na proces o zmiennej produktywności pozostawiamy jako temat przyszłych badań.
4.1.4 Model ocen i.i.d. Przegląd rozkładów
Przedyskutujemy teraz model ocen i.i.d., w którym najłatwiej wyprowadza się wyniki teoretyczne. Rozpatrujemy w nim losowe oceny ustalonych z góry n > 0 produktów. Oceny te są reprezentowane przez niezależne zmienne losowe (X1, . . . , Xn) pochodzące z pewnego rozkładu F określonego na prze-dziale I = [a, b].
Interesuje nas szczególnie przypadek, gdy F jest ciągły, gdyż może to Rozkłady ciągłe
znacząco uprościć analizę. Jak zobaczymy dalej, założenie to często jest ade-kwatne, nawet gdy w danym obszarze zainteresowań oceny produktów po-chodzą ze zbioru dyskretnego. W wielu zastosowaniach przyjmujemy, że gę-stość f rozpatrywanego rozkładu spełnia warunki f′(x) < 0 oraz f′′(x) > 0, a więc jest to, odpowiednio, funkcja malejąca i wypukła (por. omówione wy-żej prawo Lotki). Takie założenie wynika z natury niektórych problemów, w których wyższa ocena jest trudniejsza do zdobycia niż niższa, przy czym prawdopodobieństwo uzyskania wysokiej oceny maleje szybciej niż liniowo.
Niektóre tego typu rozkłady mogą się nawet cechować tzw. ciężkimi ogonami.
Przykładami takich rozkładów, określonych na I = R+, jest rozkład wy-kładniczy bądź rozkład Pareto.
Analogiczne warunki można rozpatrywać dla funkcji masy prawdopodo-bieństwa w przypadku rozkładów dyskretnych.
W niniejszej pracy nie będziemy rozważać przypadku zmiennych losowych Niezależność
nie spełniających warunku niezależności. Taki model mógłby być przydatny, gdybyśmy brali na przykład pod uwagę zwiększającą się w czasie „renomę”
producentów, która mogłaby powodować dodatkowy wzrost ocen kolejnych wytwarzanych produktów. To interesujące lecz skomplikowane zagadnienie (problemem samym w sobie jest już dobór postaci zależności między zmien-nymi) pozostawiamy jako kolejny temat do przyszłych badań.
Poniżej przyjrzymy się dwóm rodzinom rozkładów najczęściej polecanym w literaturze do modelowania interesującej nas klasy zjawisk, w tym biblio-metrycznych. Będą to: rozkład Pareto II rodzaju oraz jego zdyskretyzowana wersja.
Rozkład Pareto II rodzaju P2(k, s)
W literaturze można spotkać wiele uogólnień klasycznego rozkładu Pareto (GPD, ang. Generalized Pareto Distributions), zob. np. [184, 198]. Tutaj omówimy tzw. rozkład Pareto II rodzaju, zwany czasem rozkładem Lomaxa.
Powiemy, że zmienna losowa X pochodzi z rozkładu Pareto II rodzaju o parametrach kształtu k > 0 i skali s > 0, co oznaczamy przez X ∼ P2(k, s), jeżeli jej rozkład jest opisany gęstością:
Rodzina P2(k, s)
f (x) = ksk
(s + x)k+1 (x 0). (4.4)
Dystrybuanta tego rozkładu dana jest wzorem:
F (x) = 1− sk
(s + x)k (x 0). (4.5)
Istotną własnością parametru kształtu jest to, że jeśli X ∼ P2(kx, s) oraz Y ∼ P2(ky, s), gdzie kx < ky, to X dominuje Y stochastycznie, co oznaczamy przez X st Y.
Relacja st
Z drugiej strony, jeśli zmienne losowe pochodzą z rozkładów o tych samych parametrach kształtu, to większy parametr skali implikuje stochastyczną dominację zmiennej o mniejszej wartości s, tzn. jeśli X ∼ P2(k, sx) oraz Y ∼ P2(k, sy), to dla sx > sy zachodzi X st Y.
Wartość oczekiwana zmiennej losowej X ∼ P2(k, s) jest określona dla k > 1 i równa jest
EX = s k− 1. Wariancja istnieje dla k > 2 i wynosi
Var X = ks2
(k− 2)(k − 1)2. Ogólnie, i-ty moment zwykły jest skończony dla k > i:
EXi = Γ(i + 1)Γ(k− i) Γ(k + 1) ksi.
W bibliotece CITAN udostępniliśmy funkcje związane z omawianym roz-kładem:dpareto2()podaje wartość gęstości (4.4),ppareto2()zwraca war-tość dystrybuanty (4.5), a qpareto2() oblicza funkcję kwantylową (por.
rozdz. B.4.4).
Generowanie liczb losowych. Do utworzenia generatora liczb losowych z rozkładu Pareto II rodzaju możemy wykorzystać metodę odwracania dys-trybuanty, gdyż funkcję kwantylową daje się wyznaczyć w sposób analityczny:
F−1(y) = s(1− y)−1/k− 1.
Generator ten został zaimplementowany w funkcji rpareto2() (por.
rozdz. B.4.4).
Estymacja punktowa parametrów. Estymatory parametrów rozkładu Pareto II rodzaju można próbować wyznaczyć metodą największej wiarogod-ności (MLE). Logarytm funkcji wiarogodwiarogod-ności dany jest wzorem
lnL(k, s; x) = n ln k + nk ln s − (k + 1)
Dla n > 2 zachodzi
Ebk = k n n− 1, Varbk = k2 n2
(n− 2)(n − 1)2,
jest to zatem tylko asymptotycznie nieobciążony estymator parametru k.
Łatwo jednak uwzględnić poprawkę na nieobciążoność: Nieobciążony estymator MLE:
Wartość tego estymatora można wyznaczyć za pomocą funkcji pareto2.mle-kestimate(), zob. rozdz. B.4.4.
W przypadku nieznanego parametru skali, zakładając, że lokalne maksi-mum funkcji wiarogodności istnieje (co w wielu przypadkach nie musi być prawdą, zob. szczegółową analizę tego zagadnienia w [54]), s i k wyznaczamy z układu równań Zauważmy, że drugie równanie jest zależne tylko od zmiennej s, może być zatem rozwiązane np. za pomocą jednowymiarowej metody Newtona-Raphsona. Potem wynikowe sb wystarczy podstawić do pierwszego równania
w celu otrzymania wartości k. Estymatory te, oznaczone dalej przezb bkMLE Estymatory
bkMLE ibsMLE
i sbMLE, można wyznaczyć za pomocą funkcji pareto2.mleksestimate(), zob. rozdz. B.4.4.
W przypadku niemożliwości obliczenia wartości estymatorów metodą naj-większej wiarogodności powinniśmy użyć innego estymatora. Metoda momen-tów nie jest polecana, gdyż momenty zwykłe rzędu i, jak wspomniano wyżej, są określone tylko dla k > i. Co więcej, w przypadku niektórych zaobserwo-wanych wartości x, również estymatory wyznaczone metodą kwantyli mogą nie być wyznaczalne. Inne, bardziej skomplikowane metody zostały opisane m.in. w pracach [17, 126, 198]. Jak pokazują wyniki symulacyjne, dla nie-znanego zarówno parametru k jak i s, estymatory te cechują się lepszymi własnościami niż estymator największej wiarogodności.
Warto zwrócić szczególną uwagę na metodę szacowania wartości parame-trów zaproponowaną przez Zhanga i Stevensa [199]. Przedstawiany estymator cechuje się niewielkim obciążeniem (najczęściej dodatnim), względnie małą wariancją oraz tym, że da się wyznaczyć jego wartość w prawie wszystkich przypadkach.
Niech θ := −1s, gdzie −1 ¬ θ < 0. Poszukujemy estymatora bay-esowskiego minimalizującego średniokwadratową funkcję ryzyka (estymator MMSE, ang. minimum mean square error estimator), który jest wartością oczekiwaną rozkładu a posteriori dla pewnej gęstości rozkładu a priori h(θ; x) oraz logarytmu z profilowej funkcji wiarogodności (ang. profile log-likelihood) parametru θ postaci
ℓ(θ; x) := n (ln(−k∗θ)− 1 − 1/k∗) ,
W rozpatrywanej pracy stosowane jest następujące przybliżenie nume-ryczne wzoru (4.10). Dla m = 20 + ⌊n⌋ i j = 1, . . . , m niech
Ponadto, niech w(ϑj) = exp(ℓ(ϑj; x))/Pmi=1exp(ℓ(ϑi; x)). W konsekwencji, estymatorem parametru skali będzie sbZS = −1/Pmj=1ϑjw(ϑj). Wartość kbZS
wyznaczamy z pierwszego równania układu (4.9). Jako rozkładu a priori Estymatory
MMSE:
bkZSibsZS
h(θ; x) używamy tutaj gęstości rozkładu P2(2, 1/3 x(⌊n/4+0.5⌋)) (por. mo-tywację i dyskusję w [199]). Algorytm wyznaczający estymatory przed-stawioną metodą został zaimplementowany w bibliotece CITAN w funkcji pareto2.zsestimate(), zob. rozdz. B.4.4.
Przykład 4.1.1. Rys. 4.2–4.3 przedstawiają obciążenie i odchylenie stan-dardowe estymatorów kbMLEU (dla znanego s), bkMLE, bkZS, sbMLE i sbZS dla roz-kładów P2(1, 1), P2(5, 0,5) oraz P2(2, 2) uzyskane za pomocą M = 10000 symulacji Monte Carlo. Niestety, widzimy, że brak informacji o parametrze skali sprawia, że błąd oszacowania znacząco wzrasta, a dla k ≫ 1 estymatory Zhanga-Stevensa (a tym bardziej największej wiarogodności) zachowują się
wręcz słabo.
Testowanie zgodności. Przedstawione estymatory mogą być użyte w pro- GoF cedurze testowania zgodności danej zaobserwowanej próbki losowej z roz-kładem Pareto II rodzaju. W tym celu można użyć np. testu Kołmogo-rowa, Craméra-von Misesa bądź Andersona-Darlinga (który przykłada więk-szą wagę od pozostałych do obserwacji w ogonach rozkładu), por. analizę własności tych testów w pracach [46, 199].
W bibliotece CITAN udostępniliśmy funkcję funkcję pareto2.goftest() (zob. rozdz. B.4.4), która domyślnie stosuje test Andersona-Darlinga, w któ-rym rozkład statystyki przybliżony jest za pomocą metody zaproponowanej w [137]. To postępowanie (wraz z omówioną wyżej metodą estymacji) jest równoważne z przedstawionym w pracy [199].
Zastosowanie. Rozkład Pareto II rodzaju stosuje się często do modelo-wania tzw. wartości ekstremalnych [por. 46]. Jest on zalecany w literaturze bibliometrycznej m.in. do opisu rozkładu liczby cytowań, por. np. [93] oraz prace [11, 91, 94]. Często zakłada się ponadto, iż s 1.
Interesującą własnością tego rozkładu jest to, że jeśli X ∼ P2(k, s), to rozkład warunkowy zmiennej X −t pod warunkiem X > t, dla każdego t 0, jest rozkładem P2(k, s + t).
20 40 60 80 100
0.00.20.40.60.81.01.2
P2(1, 1); M=10000
n
Obciazenie
k_MLEU k_MLE k_ZS s_MLE s_ZS
20 40 60 80 100
0.00.51.01.52.02.53.0
n
Odch.std.
k_MLEU k_MLE k_ZS s_MLE s_ZS
Rysunek 4.2: Obciążenie i odchylenie standardowe różnych estymatorów pa-rametrów rozkładu P2(1, 1).
20 40 60 80 100
0.00.51.0
P2(2, 2); M=10000
n
Obciazenie
k_MLEU k_MLE k_ZS s_MLE s_ZS
20 40 60 80 100
0.51.01.52.02.53.03.5
n
Odch.std.
k_MLEU k_MLE k_ZS s_MLE s_ZS
Rysunek 4.3: Obciążenie i odchylenie standardowe różnych estymatorów pa-rametrów rozkładu P2(2, 2).
20 40 60 80 100
−2−1012345
P2(5, 0.5); M=10000
n
Obciazenie
k_MLEU k_MLE k_ZS s_MLE s_ZS
20 40 60 80 100
051015
n
Odch.std.
k_MLEU k_MLE k_ZS s_MLE s_ZS
Rysunek 4.4: Obciążenie i odchylenie standardowe różnych estymatorów pa-rametrów rozkładu P2(5, 0,5).
Zdyskretyzowany rozkład Pareto II rodzaju DP2(k, s)
W pewnych zastosowaniach rozpatruje się naturalną dyskretyzację rozkładu Rodzina
DP2(k, s)
Pareto II rodzaju, tj. zmienną losową ⌊X⌋ ∼ DP2(k, s), gdzie X ∼ P2(k, s).
W tym przypadku funkcja masy prawdopodobieństwa dla i = 0, 1, 2, . . . jest postaci
P(⌊X⌋ = i) = sk(s + i)−k− (s + i + 1)−k. (4.12) Przypomnijmy, że indeks h można uzyskać za pomocą przekształcenia funkcją ⌊·⌋ pewnej S-statystyki, por. tw. 2.4.23.
Funkcjediscrpareto2.mlekestimate()oraz discrpareto2.mlekses-timate() (zob. rozdz. B.4.4) implementują procedury wyznaczające nume-ryczne przybliżenie wartości estymatorów największej wiarogodności. Po-nadto,discrpareto2.goftest() implementuje test zgodności chi-kwadrat.
Przykład 4.1.2. Rys. 4.5 i 4.6 pokazują podstawowe własności omówionych estymatorów parametrów rozkładu DP2(2, 2) (zwróćmy uwagę na skalę lo-garytmiczną na osi y) oraz DP2(1,5, 1,5) (duża liczność prób) uzyskane za pomocą M = 5000 iteracji MC. Zauważmy, że dla niewielkich n w przypadku nieznanego s estymator największej wiarogodności sprawuje się bardzo słabo.
W dalszej części niniejszego rozdziału większość przykładów będziemy
ilustrować za pomocą ciągłych rozkładów z rodziny Pareto II rodzaju (w mo- Wybór
rozkładów do ilustracji
delu ocen i.i.d.) lub ewentualnie jego zdyskretyzowanej wersji. Jak wspo-mnieliśmy, założenie ciągłości upraszcza obliczenia i wyprowadzanie wyni-ków teoretycznych. Ponadto, jak wskazują wyniki z rozdziału 6, stosowanie tej rodziny jest także uzasadnione empirycznie w najbardziej interesującym nas obszarze zastosowań, tj. bibliometrii.
20 40 60 80 100
0.020.100.502.0010.00
DP2(2, 2); M=5000
n
Obciazenie
k_MLEk k_MLEks s_MLEks
20 40 60 80 100
0.20.51.02.05.0
n
Odch.std.
k_MLEk k_MLEks s_MLEks
Rysunek 4.5: Obciążenie i odchylenie standardowe wybranych estymatorów parametrów rozkładu DP2(2, 2).
2000 4000 6000 8000 10000
0.0000.0100.0200.030
DP2(1.5, 1.5); M=1000
n
Obciazenie
k_MLEk k_MLEks s_MLEks
2000 4000 6000 8000 10000
0.020.050.100.20
n
Odch.std.
k_MLEk k_MLEks s_MLEks
Rysunek 4.6: Obciążenie i odchylenie standardowe wybranych estymatorów parametrów rozkładu DP2(1,5, 1,5).