Besoin d'aide programmation vraisemblance

Bonjour,

Je vous sollicite aujourd'hui parce que j'ai un souci avec ma vraisemblance.
Sans rentrer dans les détails, le problème est que mes paramètres ne peuvent pas prendre certaines valeurs: j'ai des logs, et des sqrt (racines carrés) qui me restreignent certains paramètres à être positifs par ex.

Donc ma question est: comment faites-vous dans ce cas pour signifier à R de ne pas prendre ces valeurs lorsque l'algorithme d'optimisation est en route?

J'utilise la fonction maxLik et les algo BFGS et NM.

Merci d'avance.

Il serait peut-être utile que tu mettes ici ce que tu as déjà fait, pour mieux voir là où ça bloque

Bonjour,

Pour augmenter vos chances d'obtenir de l'aide, il est toujours bienvenu de donner un exemple de code reproductible avec, lorsque le cas s'y prête, un (court) exemple de données sur lequel appliquer le code et les éventuels messages d'erreurs.

De plus, lorsque vous citez les fonctions R utilisées (par exemple "maxLik" ici), il est préférable aussi de préciser le package auquel elles appartiennent.

Bonne continuation :)

Cordialement,

A.D.

1 pièce(s) jointe(s)

Oui vous avez raison. La base de données est en pièce jointe.
Voici le code en question:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
data_AP <- read.table("/data_AP.txt",header=T,sep="\t",quote="")
ap<-data_AP
names(ap)
summary(ap)
attach(ap)
 
age[is.na(age)]<-0
concern[is.na(concern)]<-0
y1<-(yes1=="1")*1
sex<-(sex=='1')*1
inc1000<-income/1000
 
library(maxLik)
library(MASS)
 
 
 
#Reg 4: SML, DCCV estimation
######Here sigma is ,ormalized to 1, pa[2] is theta, not sigma
 
 
approbit<- glm(y1 ~ educ + age + sex + concern, family=binomial(link="probit"), data=ap)
summary(approbit)
 
set.seed(123)
tau<-rnorm(500)
 
 
L4<-function(par){
	pa<-par[1:8]
 
## Là je contrains mes paramètres pour que la vraisemblance puisse être calculée: 
## racine carrée et log doivent être sur des nombres positifs.
## pa[8][pa[8]<0]<-NaN  veut dire que le paramètre 8 devient NaN quand il est négatif,
## du coup sqrt(pa[8]) peut maintenant être calculée.
 
 
	pa[8][pa[8]<0]<-NaN
	pa[8][pa[8]>=1]<-NaN
 
	pa[1][pa[1]>=-pa[8]]<-NaN
 
	pa[1][pa[1]<=-1]<-NaN
	pa[1][pa[1]>=0]<-NaN
 
 
	u0<-tau
 
	alpha2s<-pa[3]*educ+pa[4]*age+pa[5]*sex+pa[6]*concern
 
	om<- list()
	for (i in 1:500)
	{
	om[[i]]<-c(exp(alpha2s-u0[i])
	+exp(pa[2])*(travcost-pt1)
	+exp(pa[7])*pt1)
 
	om[[i]][om[[i]]<=1e-200]<-1e-200
	}
 
 
 
	omega<- list()
	for (i in 1:500)
	{
	omega[[i]]<-c(log(om[[i]])
	+(1+pa[1])*(log(1+pa[1])-pa[2]-log(travcost))
	+pa[1]*alpha2s-pa[1]*log(-pa[1]))
	}
 
	OMEGA<-do.call(cbind,omega)
 
 
	SP<-rowMeans(pnorm((OMEGA-pa[8]*u0)
	/sqrt(pa[1]^2 - pa[8]^2)))
	mSP<-1-SP
 
	SP[SP<=1e-200]<-1e-200
	SP[SP>=1]<-1
 
    mSP[mSP<=1e-200]<-1e-200
	mSP[mSP>=1]<-1
 
	S<-sum(y1*log(SP)+(1-y1)*log(mSP))
 
 
	return(S)
	}
 
start4<-c(-0.1,0.1,0.069,-0.015,0.372,0.468,0.1,0.01)
 
fstart4_100sim<-c(-0.7028468, -1.046796, 0.2159982, -0.03577681, 0.7481096, 1.185199, -13.85324, 0.3909023)
fstart4_500sim<-c(-0.6954158, -1.049007, 0.207197, -0.03550529, 0.7908098, 1.142287, -13.74835, 0.4854157)
 
st4<-rnorm(8)
 
opt4<-maxLik(L4,start=start4, print.level=2,method='NM')
 
opt4
summary(opt4)
 
sta4<-c(coef(opt4))