[R] optimize simultaneously two binomials inequalities using nlm( ) or optim( )

Marc Schwartz marc_schwartz at comcast.net
Tue Aug 5 18:11:35 CEST 2008


Just a quick follow up to Spencer's post, you might want to look at the 
AcceptanceSampling package on CRAN:

http://cran.r-project.org/web/packages/AcceptanceSampling/index.html

HTH,

Marc Schwartz


on 08/05/2008 09:00 AM Spencer Graves wrote:
>      I saw your post on 7/29, and I have not seen a reply, so I will 
> attempt a response to the question at the start of your email: 
> obtain the smallest value of 'n' (sample size)
> satisfying both inequalities:
>             (1-alpha) <= pbinom(c, n, p1) && pbinom(c, n, p2) <= beta,
> where alpha, p1, p2, and beta are given, and I assume that 'c' is also 
> given, though that's not perfectly clear to me. 
>       Since 'n' is an integer, standard optimizers like optim and nlm 
> are not really appropriate.  This sounds like integer programming.  
> RSiteSearch('integer programming', 'fun') just produced 138 hits for 
> me.  You might find something useful there. 
>       However, before I tried that, I'd try simpler things first.  
> Consider for example the following: 
> c. <- 5
> # I used 'c.' not 'c', because 'c' is the name of a function in R.  
> alpha <- .05
> beta <- .8
> p1 <- .2
> p2 <- .9
> n <- c.:20
> p.1 <- pbinom(c., n, p1)
> p.2 <- pbinom(c., n, p2)
> 
> good <- (((1-alpha) <= p.1) & (p.2 <= beta))
> min(n[good])
> 
> op <- par(mfrow=c(2, 1)) plot(n, p.1)
> abline(h=1-alpha)
> plot(n, p.2)
> abline(h=beta)
> par(op)
> 
>       In this case, n = 6 satisfies both inequalities, though n = 15 
> does not.  If min(n[good]) = Inf either no solution exists or you need 
> to increase the upper bound of your search range for 'n'. 
>       If you'd like more help, PLEASE do read the posting guide 
> 'http://www.R-project.org/posting-guide.html' and provide commented, 
> minimal, self-contained, reproducible code.
> 
>       Hope this helps.        Spencer Graves
> 
> emir.toktar at gmail.com wrote:
>> Dear R users,
>>
>> I´m trying to optimize simultaneously two binomials inequalities (used in
>> acceptance sampling) which are nonlinear solution, so there is no simple
>> direct solution. Please, let me explain shortly the the problem and the
>> question as following.
>>
>> The objective is to obtain the smallest value of 'n' (sample size)
>> satisfying both inequalities:
>> (1-alpha) <= pbinom(c, n, p1) && pbinom(c, n, p2) <= beta
>>
>> Where p1 (AQL) and p2 (LTPD) are probabilities parameters (Consumer and
>> Producer), the alpha and beta are conumer and producer risks, and 
>> finally,
>> the 'n' represents the sample size and the 'c' an acceptance number or
>> maximum number of defects (nonconforming) in sample size.
>>
>> Considering that the 'n' and 'c' values are integer variables, it is
>> commonly not possible to derive an OC curve including the both 
>> (p1,1-alpha)
>> and (p2,beta) points. Some adjacency compromise is commonly required,
>> achieved by searching a more precise OC curve with respect to one of the
>> points.
>> I´m using Mathematica 6 but it is a Trial, so I would like use R 
>> intead (or
>> better, I need it)!
>>
>> To exemplify,  In Mathematica I call the function using NMinimize passing
>> the restriction and parameters:
>>
>> /* function name "findOpt" and parameters... */
>>
>> restriction = (1 - alpha) <= CDF[BinomialDistribution[sample_n, p1], 
>> c]         &&  betha >= CDF[BinomialDistribution[sample_n, p2], c] 
>>         &&   0 < alpha < alphamax && 0 < betha < bethamax           
>> &&   1 < sample_n <= lot_Size   &&   0 <= c < lot_size
>>         &&  p1 < p2 < p2max ;
>>
>> fcost =     sample_n/lot_Size;
>> result = NMinimize[{fcost, restriction}, {sample_n, c, alpha, betha, 
>> p2max},
>> Method -> "NelderMead", AccuracyGoal -> 10];
>>
>> /* Calling the function findOpt */
>> findOpt[p1=0.005, lot_size=1000, alphamax=0.05, bethamax =0.05, p2max =
>> 0.04]
>>
>> /* and I got the return of values of; minimal "n", "c", "alpha", 
>> "betha" and
>> the "p2" or (LTPD) were computed */ {0.514573, {alpha$74 -> 0.0218683,
>> sample_n$74 -> 155.231, betha$74 -> 0.05,
>> c$74 -> 2, p2$74 -> 0.04}}
>>
>> Now, using R, I would define the "pbinom(c, n, prob)" given only the 
>> minimum
>> and maximum values to "n" and "c" and limits to p1 and p2 probabilities
>> (Consumer and Producer), similar on the example above. 
>> I found some examples to optimize equations in R and some tips, but I 
>> not be
>> able to define the sintaxe to use with that functions. Among the 
>> functions
>> that could be used to resolve the problem presented, I found the function
>> optim() that it is used for unconstrained optimization and the nlm() 
>> which
>> is used for solving nonlinear unconstrained minimization problems. May 
>> I wrong, but the nlm() function would be appropriate to solve this
>> problem, is it right?
>>
>> Can I get a pointer to solve this problem using the nlm() function or 
>> where
>> could I get some tips/example to help me, please?
>>
>> // (1-alpha) <= pbinom(c, n, p1) && pbinom(c, n, p2) <= beta
>> It was used "betha" parameter name to avoid the 'beta' function used in
>> Mathematica...
>>
>>
>> findS <- function(p1='numeric', lot_size='numeric', alphamax='numeric',
>> bethamax ='numeric', p2max ='numeric')
>> {
>>     (1 - alpha) <= pbinom(c, sample_n, p1) &&  betha >= pbinom(c,
>> sample_n, p2)
>>     &&   0 < alpha < alphamax && 0 < betha < bethamax       &&   1 < 
>> sample_n <= lot_Size   &&   0 <= c < lot_size
>>     &&  p1 < p2 < p2max ;
>> }
>>
>> Parameters:
>>     p1=0.005,     lot_size=1000,     alphamax=0.05,     bethamax 
>> =0.05,     p2max = 0.04
>>
>>
>> Minimize results should return/printing the following values:
>>     sample_n,     (minimal sample size)
>>     c ,         (critical level of defectives)
>>     alpha ,     (producer's risk)
>>     betha ,     (consumer's risk)
>>     p2max    (consumer's probability p2)
>>
>>
>> Could one help me understand how can desing the optimize nonlinear 
>> function
>> using R for two binomials or point me some tips?
>>
>>
>> Thanks in advance.
>>
>> EToktar
>>



More information about the R-help mailing list