是否有人尝试在R中复制斯蒂格勒饮食问题?
到目前为止，这是我所拥有的:
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, Dried", "Raisins, Dried", "Peas, Dried", "Lima Beans, Dried",
"Navy Beans, Dried", "Coffee", "Tea", "Cocoa", "Chocolate",
"Sugar", "Corn Syrup", "Molasses", "Strawberry Preserves"
),
unit=c(
"10 lb.", "1 lb.", "28 oz.", "8 oz.", "1 lb.", "24 oz.",
"1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.",
"1 qt.", "14.5 oz.", "1 lb.", "1 lb.", "1 doz.", "1 lb.",
"1/2 pt.", "1 lb.", "1/2 pt.", "1 lb.", "1 lb.", "1 lb.",
"1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.",
"1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.", "1 lb.",
"16 oz.", "1 lb.", "1 lb.", "1 doz.", "1 doz.", "1 lb.",
"1 lb.", "1 bunch", "1 stalk", "1 head", "1 lb.", "15 lb.",
"1 lb.", "1 lb.", "No. 2 1/2", "No. 2 1/2", "No. 2 1/2",
"No. 2", "No. 2", "16 oz.", "No. 2", "No. 2", "No. 2", "10 1/2 oz.",
"1 lb.", "1 lb.", "15 oz.", "1 lb.", "1 lb.", "1 lb.", "1 lb.",
"1/4 lb.", "8 oz.", "8 oz.", "10 lb.", "24 oz.", "18 oz.",
"1 lb."
),
price=c(
36, 14.1, 24.2, 7.1, 4.6, 8.5, 7.5, 7.1, 7.9, 9.1,
9.1, 24.8, 15.1, 11, 6.7, 30.8, 16.1, 32.6, 24.2, 14.1, 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, 24.2, 25.6, 27.4, 16, 30.3, 42.3, 13, 4.4, 6.1,
26, 30.9, 7.1, 3.7, 4.7, 7.3, 8.2, 3.6, 34, 8.1, 5.1, 16.8,
20.4, 21.3, 27.7, 10, 7.1, 10.4, 13.8, 8.6, 7.6, 15.7, 9,
9.4, 7.9, 8.9, 5.9, 22.4, 17.4, 8.6, 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, 6.1, 8.4, 10.8, 20.6, 2.9, 7.4, 3.5, 15.7,
8.6, 20.1, 41.7, 2.9, 2.2, 3.4, 3.6, 8.5, 2.2, 3.1, 3.3,
3.5, 4.4, 10.4, 6.7, 18.8, 1.8, 1.7, 5.8, 5.8, 4.9, 1, 2.2,
2.4, 2.6, 2.7, 0.9, 0.4, 5.8, 14.3, 1.1, 9.6, 3.7, 3, 2.4,
0.4, 1, 7.5, 5.2, 2.3, 1.3, 1.6, 8.5, 12.8, 13.5, 20, 17.4,
26.9, 0, 0, 8.7, 8, 34.9, 14.7, 9, 6.4
),
protein=c(
1411, 418, 377, 252,
897, 680, 460, 907, 488, 484, 439, 130, 288, 310, 422, 9,
17, 238, 448, 49, 661, 18, 0, 0, 166, 214, 213, 309, 404,
333, 245, 140, 196, 249, 152, 212, 164, 184, 156, 705, 27,
60, 21, 40, 138, 125, 73, 51, 27, 166, 336, 106, 138, 20,
8, 16, 33, 54, 364, 136, 136, 63, 71, 87, 99, 104, 1367,
1055, 1691, 0, 0, 237, 77, 0, 0, 0, 11
),
calcium=c(
2, 0.7, 14.4, 0.1,
1.7, 0.8, 0.6, 5.1, 2.5, 2.7, 1.1, 0.4, 0.5, 10.5, 15.1,
0.2, 0.6, 1, 16.4, 1.7, 1, 0.2, 0, 0, 0.1, 0.1, 0.1, 0.2,
0.2, 0.2, 0.1, 0.1, 0.2, 0.3, 0.2, 0.2, 0.1, 0.1, 0.1, 6.8,
0.5, 0.4, 0.5, 1.1, 3.7, 4, 2.8, 3, 1.1, 3.8, 1.8, 0, 2.7,
0.4, 0.3, 0.4, 0.3, 2, 4, 0.2, 0.6, 0.7, 0.6, 1.7, 2.5, 2.5,
4.2, 3.7, 11.4, 0, 0, 3, 1.3, 0, 0.5, 10.3, 0.4
),
iron=c(
365, 54, 175, 56, 99, 80, 41, 341, 115, 125, 82, 31, 50, 18, 9, 3,
6, 52, 19, 3, 48, 8, 0, 0, 34, 32, 33, 46, 62, 139, 20, 15,
30, 37, 23, 31, 26, 30, 24, 45, 36, 30, 14, 18, 80, 36, 43,
23, 22, 59, 118, 138, 54, 10, 8, 8, 12, 65, 134, 16, 45,
38, 43, 173, 154, 136, 345, 459, 792, 0, 0, 72, 39, 0, 74, 244, 7
),
vit_a=c(
0, 0, 0, 0, 30.9, 0, 0, 0, 0, 0, 0, 18.9, 0, 16.8,
26, 44.2, 55.8, 18.6, 28.1, 16.9, 0, 2.7, 0, 0.2, 0.2, 0.4,
0, 0.4, 0, 169.2, 0, 0, 0, 0, 0, 0, 0, 0.1, 0, 3.5, 7.3,
17.4, 0, 11.1, 69, 7.2, 188.5, 0.9, 112.4, 16.6, 6.7, 918.4,
290.7, 21.5, 0.8, 2, 16.3, 53.9, 3.5, 12, 34.9, 53.2, 57.9,
86.8, 85.7, 4.5, 2.9, 5.1, 0, 0, 0, 0, 0, 0, 0, 0, 0.2
),
vit_b1=c(
55.4, 3.2, 14.4, 13.5, 17.4, 10.6, 2, 37.1, 13.8, 13.9, 9.9, 2.8,
0, 4, 3, 0, 0.2, 2.8, 0.8, 0.6, 9.6, 0.4, 0, 0, 2.1, 2.5,
0, 1, 0.9, 6.4, 2.8, 1.7, 17.4, 18.2, 1.8, 9.9, 1.4, 0.9,
1.4, 1, 3.6, 2.5, 0.5, 3.6, 4.3, 9, 6.1, 1.4, 1.8, 4.7, 29.4,
5.7, 8.4, 0.5, 0.8, 2.8, 1.4, 1.6, 8.3, 1.6, 4.9, 3.4, 3.5,
1.2, 3.9, 6.3, 28.7, 26.9, 38.4, 4, 0, 2, 0.9, 0, 0, 1.9, 0.2
),
vit_b2=c(
33.3, 1.9, 8.8, 2.3, 7.9, 1.6, 4.8, 8.9, 8.5, 6.4, 3,
3, 0, 16, 23.5, 0.2, 0, 6.5, 10.3, 2.5, 8.1, 0.5, 0, 0.5,
2.9, 2.4, 2, 4, 0, 50.8, 3.9, 2.7, 2.7, 3.6, 1.8, 3.3, 1.8,
1.8, 2.4, 4.9, 2.7, 3.5, 0, 1.3, 5.8, 4.5, 4.3, 1.4, 3.4,
5.9, 7.1, 13.8, 5.4, 1, 0.8, 0.8, 2.1, 4.3, 7.7, 2.7, 2.5,
2.5, 2.4, 4.3, 4.3, 1.4, 18.4, 38.2, 24.6, 5.1, 2.3, 11.9,
3.4, 0, 0, 7.5, 0.4
),
vit_b3=c(
441, 68, 114, 68, 106, 110, 60, 64,
126, 160, 66, 17, 0, 7, 11, 2, 0, 1, 4, 0, 471, 0, 0, 5,
69, 87, 0, 120, 0, 316, 86, 54, 60, 79, 71, 50, 0, 68, 57,
209, 5, 28, 4, 10, 37, 26, 89, 9, 11, 21, 198, 33, 83, 31,
5, 7, 17, 32, 56, 42, 37, 36, 67, 55, 65, 24, 162, 93, 217,
50, 42, 40, 14, 0, 5, 146, 3
),
vit_c=c(
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 177, 60, 0, 0, 0, 0,
17, 0, 0, 0, 0, 0, 0, 0,
0, 0, 525, 0, 0, 0, 0, 0, 0, 0, 46, 0, 0, 544, 498, 952,
1998, 862, 5369, 608, 313, 449, 1184, 2522, 2755, 1912, 196,
81, 399, 272, 431, 0, 218, 370, 1253, 862, 57, 257, 136,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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]
我认为如果营养数据是按单位而不是按美元计算的话，这个问题会更有指导意义。