问题描述
有人试过在 R 中复制 The Stigler Diet problen 吗？
到目前为止，这是我所拥有的：
library(lpSolve)
library(linprog)
# where d is the nutrients data frame found in the above link
f.obj_ <- d$`1939 price (cents) `
f.con_ <- matrix(c( d$Calories,d$`Protein (g)`,d$`Calcium (g)`,d$`Iron (mg)`,d$`Vitamin A (IU)`,d$`Thiamine (mg)`,d$`Riboflavin (mg)`,d$`Niacin (mg)`,d$`Ascorbic Acid (mg)`),nrow = 9,byrow = TRUE)
f.dir <- rep(">=",9)
f.rhs <- c(3.0,70.0,0.8,12.0,5.0,1.8,2.7,18.0,75.0)
lp("min",f.obj_,f.con_,f.dir,f.rhs)
# returns: 0.7164608
我的猜测是 Python 代码使用的优化器与 lpSolve（Glop 方法）中的优化器不同。 R 中有 Glop 优化器可以让我从 Py 代码中复制结果吗？
解决方法
虽然我还没有复制这个，但我已经计划了很长时间。现在你的问题是一个很好的理由，让我远离了 Netflix。
这是使用 lpSolveAPI 的解决方案：
library(lpSolveAPI)
# minimum requirements
requirements <- c(
calories=3,protein=70,calcium=0.8,iron=12,vit_a=5,vit_b1=1.8,vit_b2=2.7,vit_b3=18,vit_c=75
)
# nutrition data
dat <- data.frame(
commodity=c(
"Wheat Flour (Enriched)","Macaroni","Wheat Cereal (Enriched)","Corn Flakes","Corn Meal","Hominy Grits","Rice","Rolled Oats","White Bread (Enriched)","Whole Wheat Bread","Rye Bread","Pound Cake","Soda Crackers","Milk","Evaporated Milk (can)","Butter","Oleomargarine","Eggs","Cheese (Cheddar)","Cream","Peanut Butter","Mayonnaise","Crisco","Lard","Sirloin Steak","Round Steak","Rib Roast","Chuck Roast","Plate","Liver (Beef)","Leg of Lamb","Lamb Chops (Rib)","Pork Chops","Pork Loin Roast","Bacon","Ham,smoked","Salt Pork","Roasting Chicken","Veal Cutlets","Salmon,Pink (can)","Apples","Bananas","Lemons","Oranges","Green Beans","Cabbage","Carrots","Celery","Lettuce","Onions","Potatoes","Spinach","Sweet Potatoes","Peaches (can)","Pears (can)","Pineapple (can)","Asparagus (can)","Green Beans (can)","Pork and Beans (can)","Corn (can)","Peas (can)","Tomatoes (can)","Tomato Soup (can)","Peaches,Dried","Prunes,"Raisins,"Peas,"Lima Beans,"Navy Beans,"Coffee","Tea","Cocoa","Chocolate","Sugar","Corn Syrup","Molasses","Strawberry Preserves"
),unit=c(
"10 lb.","1 lb.","28 oz.","8 oz.","24 oz.","1 qt.","14.5 oz.","1 doz.","1/2 pt.","16 oz.","1 bunch","1 stalk","1 head","15 lb.","No. 2 1/2","No. 2","10 1/2 oz.","15 oz.","1/4 lb.","10 lb.","18 oz.","1 lb."
),price=c(
36,14.1,24.2,7.1,4.6,8.5,7.5,7.9,9.1,24.8,15.1,11,6.7,30.8,16.1,32.6,17.9,16.7,20.3,9.8,39.6,36.4,29.2,22.6,14.6,26.8,27.6,36.6,30.7,25.6,27.4,16,30.3,42.3,13,4.4,6.1,26,30.9,3.7,4.7,7.3,8.2,3.6,34,8.1,5.1,16.8,20.4,21.3,27.7,10,10.4,13.8,8.6,7.6,15.7,9,9.4,8.9,5.9,22.4,17.4,16.2,51.7,13.7,13.6,20.5
),calories=c(44.7,11.6,11.8,11.4,36,28.6,21.2,25.3,15,12.2,12.4,8,12.5,8.4,10.8,20.6,2.9,7.4,3.5,20.1,41.7,2.2,3.4,3.1,3.3,18.8,1.8,1.7,5.8,4.9,1,2.4,2.6,2.7,0.9,0.4,14.3,1.1,9.6,3,5.2,2.3,1.3,1.6,12.8,13.5,20,26.9,8.7,34.9,14.7,6.4
),protein=c(
1411,418,377,252,897,680,460,907,488,484,439,130,288,310,422,17,238,448,49,661,18,166,214,213,309,404,333,245,140,196,249,152,212,164,184,156,705,27,60,21,40,138,125,73,51,336,106,33,54,364,136,63,71,87,99,104,1367,1055,1691,237,77,11
),calcium=c(
2,0.7,14.4,0.1,0.8,0.6,2.5,0.5,10.5,0.2,16.4,0.3,6.8,4,2.8,3.8,2,4.2,10.3,0.4
),iron=c(
365,175,56,80,41,341,115,82,31,50,6,52,19,48,32,46,62,139,30,37,23,24,45,14,43,22,59,118,12,65,134,38,173,154,345,459,792,72,39,74,244,7
),vit_a=c(
0,18.9,44.2,55.8,18.6,28.1,16.9,169.2,11.1,69,7.2,188.5,112.4,16.6,918.4,290.7,21.5,16.3,53.9,53.2,57.9,86.8,85.7,4.5,0.2
),vit_b1=c(
55.4,3.2,10.6,37.1,13.9,9.9,2.1,6.4,18.2,1.4,4.3,29.4,5.7,8.3,1.2,3.9,6.3,28.7,38.4,1.9,vit_b2=c(
33.3,8.8,4.8,23.5,6.5,50.8,5.4,7.7,18.4,38.2,24.6,11.9,vit_b3=c(
441,68,114,110,64,126,160,66,7,471,5,120,316,86,79,57,209,28,89,198,83,42,67,55,162,93,217,146,3
),vit_c=c(
0,177,525,544,498,952,1998,862,5369,608,313,449,1184,2522,2755,1912,81,399,272,431,218,370,1253,257,0
),stringsAsFactors=FALSE
)
rownames(dat) <- dat[,"commodity"]
dat[,"commodity"] <- NULL
# formulate problem
lprec <- lpSolveAPI::make.lp(ncol=nrow(dat))
lpSolveAPI::lp.control(lprec,sense='min')
lpSolveAPI::set.objfn(lprec,rep(1,nrow(dat)))
for (req in names(requirements)) {
lpSolveAPI::add.constraint(lprec,dat[,req],">=",requirements[req])
}
# solve
status <- lpSolveAPI::solve.lpExtPtr(lprec)
stopifnot(status==0)
message("cost per year ($): ",lpSolveAPI::get.objective(lprec)*365.25)
diet <- lpSolveAPI::get.variables(lprec)
names(diet) <- rownames(dat)
diet <- diet[diet>0]
如果营养数据是按单位而不是按美元计算，我认为这个问题会更有启发性。